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Ph- Abstract 

Basing on analogy between the three-body scattering problem and the diffraction 
| problem of the plane wave (for the case of the short range pair potentials) by the system 

i i ' of six half transparent screens, we presented a new approach to the few-body scattering 

problem. The numerical results have been obtained for the case of the short range 
nonnegative pair potentials. The presented method allows a natural generalization to 
Q\ ' the case of the long range pair potentials. 

(N : 

in ■ 

^ ; 1 Introduction 

o : 

. The quantum system of two particles interacting via the Coulomb potential is proba- 

bly the most known model of the Quantum mechanics. The model allows an explicit 
solution. Oppositely, the mathematical status of the system of three quantum particles 
with the pair Coulomb interaction is relatively poor. The system of three particles with 
short range pair interactions was successfully studied by L. Faddeev pQ, but the direct 
generalization to the Coulomb type potentials was found impossible. Something, how- 
ever, is known: the quantative nature of the spectrum and the asymptotic behavior of 
the solutions of the non-stationary Schrodinger equation. These results were obtained 
in frameworks of a non-stationary approach, see [2j [3] . Nevertheless a mathematically 
consistent stationary approach similar to the Lippmann-Schwinger integral equation, 
or something analogous, was not developed though. Such an approach is needed if 
we are interested in numerical parameters of many important physical processes like 
dissociative recombination in atomic and molecular physics with applications to as- 
trophysics, formation and break up processes of large molecules in bioengineering and 
medicine, formation of the molecular resonance states in chemical physics, dynamics 
of the few electron systems in wave-conductor nano-technology. 
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There are specific difficulties that are characteristic for the systems with Coulomb 
type interactions. 

They are naturally explained by the fact that the long range interactions crucially 
affects the asymptotic behavior at infinity in the configuration space of the eigenfunc- 
tions, Green's functions and other similar objects. The consequences of that affection 
on the structure of asymptotics up to now have not been taken in account in correct 
mathematical manner. As a result, such approaches to many particle scattering as 
Faddeev's equations AGS equations [5], successfully applicable to the systems with 
short range potentials, do not work for the systems with the long range potentials. 

The asymptotic behavior of the wave functions for the systems of few charged 
particles has been studied only in some domains of configuration space but not for 
all asymptotic directions. Let us shortly list some known results. In (6j [7] there was 
studied the asymptotic behavior of three charged particles wave function for the case 
of large distances between all three particles. Another limiting case, considered in 
[8], corresponds to configurations with one Jacobi coordinates been much larger than 
another one. 

In the list of the literature reflecting the theoretical aspects of the problem we 
mention also [TO] , [H] , [TO] , [TO] , [HI [TOJ [TO] . Application to computational aspects of 
the problem were treated in [13 [El [TO] , [H1EDI22], [23]. 

One of the typical computational approaches to such systems is to replace the 
Coulomb potentials by the Yukava potentials (or some other cut-off potentials), to 
compute the parameters of the scattering for this modified system, and to consider the 
results for small screening parameter. Mathematically, it is not a completely satisfac- 
tory procedure. Some other approximate approaches also exist. 

1.2 

We started a description of a new approach to the mathematically consistent stationary 
treatment of the scattering in the systems of three quantum particles with long range 
pair interactions in [27J. We are going to consider in turn the case of three one- 
dimensional particles with short range interaction (it is already completed theoretically 
and published), the case of three one-dimensional particles with long range (Coulomb 
type) interactions (it is also completed at the moment but is not published yet) and the 
case of three three-dimensional particles. Each next case will be based on the results 
for the preceding stage. We hope that we will be able to illustrate each theoretical 
stage by numerical computation of the field. This paper contains the numerical results 
illustrating the formulas of paper [27J. 

We assume here that the pair potentials are non- negative. In this case the spectrum 
is purely continuous, covers the positive semi-axis and is in the natural sense homo- 
geneous. In fact, this case is the most interesting at the present stage since the lower 
spectral branches for negative total energy in case of charged particles was already 
treated in |28j . 

It is worth mentioning that the scattering in the system of one-dimensional particles 
is not just a first step on the way to the case of three-dimensional particles. It is 
interesting by itself, the systems of three one-dimensional particles (neutral or charged) 
were intensively studied during many years (see, for example, j29j EUJ EU E2j [33] ) . In 
recent years there appeared a new interest to such systems since they were realized 
experimentally (see [Ml l35| l36ll3Tj ) . 

The main idea is to suggest a priori explicit formulas for the asymptotic behavior 
of the eigenfunctions of the continuous spectrum (for example, the scattered plane 
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waves) . 

The formulas describe the eigenfunctions at infinity up to the simple diverging waves 
with smooth amplitudes. If we are able to find such asymptotic behavior (satisfying 
certain criteria that will be discussed later on) even heuristically, we obtain a way 
for regular numerical computations of the eigenfunctions. We obtain simultaneously 
also a method to construct an appropriate integral equation of the same nature as 
the Lippmann-Schwinger equation for the scattering of the plane wave by a quickly 
decreasing potential that can be used to justify the asymptotic behavior rigorously 
following the ideas of [1]. 

For one-dimensional particles with quickly decreasing at infinity pair potentials 
we can use, for the description of the mentioned asymptotic behavior, the analogy 
between the stated problem and the classical problem of the diffraction of the plane 
waves by the set of semi-transparent infinite screens. This analogy was already used 
in [2i\ [25l [26l [27] . In case of long range potentials we are able to treat the diffraction 
problem analogously with the replacement of the classical plane waves by plane waves 
that are appropriately deformed by the long range tails of the Coulomb potentials. 
It is important to mention that the diffraction itself and the corresponding scattering 
problems cannot be completely reduced to the scattering of the plane waves by the 
screens; we have to add to these processes some genuine diffraction components that 
have more complicated analytical structure but still explicit description. This more 
complicated structure is also dictated by the analogy with the classical diffraction 
theory. 

Here we consider a system of three identical one-dimensional quantum particles 
interacting via short-range pair potentials. These strict limitations allow to simplify 
the narration and the view of the formulas, but the essence of the main questions 
which we are interested in and their treatments is not affected. In the following parts 
we consequentially will get rid of these limitations. As we have mentioned above, the 
theoretical part of this work is already published, but we decided for the completeness 
to repeat shortly the main theoretical ideas of (27J- The main goal of the work is to 
confirm that the approach works for the numerical computation of the eigenfunctions of 
the continuous spectrum. The approach is new even for the short range pair potentials. 

The structure of the work is as following: it consists of two parts. The first part 
is devoted to the known theoretical constructions. The second one is original and 
represents the results of numerical computer computations. 

2 Main formulas 
2.1 Configuration plane 

The configuration space of the system after the separation of motion of the center 
of mass is the hyperplane T = {x = (x\ , xi , X3) : X\ + X2 + X3 = 0} in R 3 . The 
Schroedinger equation has the form: 

- Aip + (v(x 1 )+v(x 2 )+v(x 3 ))ip = ETp, (I) 

where ip = ip(x) E C, A is the Laplace operator on T that will be described more 
specifically later on. The real- valued function v(x), x S R, is the potential of the pair 
interaction. In the present text it is supposed to be an even function with a compact 
support, v (x) = 0, |x| > b/2. We suppose E > 0. 
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The scalar product on T is given by the formula 

2 

< x,x' >= -(xixi +x 2 x' 2 + X3X3). (2) 

As usual, the norm of the vector is defined by the formula |x| 2 =< x, x > . The 
Laplacian is also generated by this scalar product. 

Let us consider on T three straight lines lj = {x : Xj = 0}, j = 1,2,3, and three 
unit vectors lj that belong to these lines and oriented such that Xj + \ increases along 
lj. Consider also the unit vectors kj that are orthogonal to lj and oriented along the 
direction of increasing of Xj. Consider, at last, three pairs of the cartesian coordinates 
(xj,yj) with respect to the bases (kj,lj). These are, so called, Jacobian coordinates 
on T. With these coordinates 

<x,x' >= xjx'j +yjy'j, \x\ 2 = x 2 + y 2 , j = 1,2,3, (3) 

and 

d 2 d 2 

A = ^ + ^, . = 1,2,3. (4) 

The lines lj define on the plane T six sectors. The internal part of a certain one 
consists of the vectors (xi,X2,X3) whose coordinates satisfy the condition Xj 1 > Xj 2 > 
Xj 3 where a = (ji, j'2, J3) is a permutation of the numbers (123). We will denote any 
sector by the corresponding permutation a and will write A = X a (see Figure 1). 




Figure 1 



Let the group S3 of permutation acts on T so that 

((7,X) -> CJX = {x h ,Xj 2 ,X j3 ), <T = (jl,32,33), X = (X!,X2,X 3 ). 
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The group contains 6 elements. The permutation can be identical, or a transposition 
of two elements, or a composition of two transpositions, some of the compositions 
coincide. Introduce the notations for the transpositions: t\ = (132), t 2 = (321), 
T3 = (213), and notice that rf = I, i = 1, 2, 3. The action of the transposition on Y 
will be denoted by the same symbol r,-, j = 1,2,3. It corresponds to the reflection 
with respect to the line lj, j = 1, 2, 3. It is clear that 

ti(x!,x 2 ,x 3 ) = (-xi,-X3,-x 2 ), n(y 1 ,y 2 ,y 3 ) = (y 1 ,y 3 ,y 2 ), (5) 

the analogous formulas are also satisfied for t 2 ,t 3 . The composition of two transpo- 
sitions generates a rotation, and the following equalities, in particular, hold: t\t 2 = 

T 3 T\ = T 2 T 3 , T 2 T\ = T1T3 = r 3 r 2 . 

Six elements r of the group S3 generate six vectors rq. If q G A CT then rq G X TU . 



2.2 Separation of variables 

Consider now the eigenfunction that describes the scattering in the system where just 
one of three potentials is not equal to zero. Now we deal with the Schrodinger equation 

- Axj +v(xj)xj = Exj. (6) 

It allows the separation of variables: 

V; U.q;. \{.rj. /.;,!( 1 (7) 

The sense of the variables (xj,yj) is clear, (kj,pj) are the Jacobian coordinates of a 
given vector q. 

The function x( x i k),x, k € R, is a solution of the ordinary differential equation 

- Xxx + v(x)x = k 2 x, (8) 

it has to be described separately. For k > there exists and is unique the solution that 
is characterized by the following asymptotic behavior: 

X(x, k) ~ s(k)e ikx , x -> +00; x(x, k) ~ e ikx + r(k)e~ ikx , x -> -00. (9) 

On the whole axis k this solution, due to the evenness of the potential, has to be 
extended by the formula x( x > k) = x{~ x i —k). Here s and r are some complex-valued 
functions of k that are called the transition and the reflection coefficients. 

We will suppose here that v(x) > therefore the equation (jHJ) does not have the 
bound states. 



2.3 Formal setting of the problem 

Our final goal is to construct the solution ip(x, q) of the Schrodinger equation that is 
characterized by the following behavior at infinity: 

e -»M|qi e »M|q| / 1 \ x 

? /; = n(x,q)^ T7 ^+/(x,q)^ T7 ^+ ^^ T ^J , x = — . (10) 

Here 
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™(x,q) = J-T-r8(x,q), (11) 
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and the 5 function has to be considered with respect to the angle measure on the unit 
circle. 

The asymptotic behavior has to be treated in a weak sense (in sense of distributions) 
with respect to x. The coefficient n before the converging circle wave coincides with the 
analogous coefficient before the converging wave in the weak asymptotic representation 
of the plane wave e* <x,q> . Therefore the solution ?/>(x, q) can be naturally called the 
scattered plane wave. 

Due to the symmetries of the potential ip(x., q) = ^(cx, crq), a G S, we always can 
assume that q belongs to a certain sector, say A/ = Ai23- We restrict ourselves here 
by the assumption that q does not belong to neighborhoods of the boundaries of the 
sector. It would be not hard to consider also the case when q belongs to the lines lj 
and their neighborhoods. 

The function / is a singular distribution. We will see that it has singularities 
on all six directions <rq, a £ S. Four of them are of 5 - function type, two (for 
crq = T2T 3q, t 2T iq) are of type of Cauchy's limiting kernel. It is worth to notice that 
although the asymptotic behavior is singular the solution itself is, naturally, a smooth 
function. 

In the case of the scattering by a quickly decreasing at infinity potential the asymp- 
totic behavior is given by the formula 

V>(x,q) = e*<^> + /(x,q)^ + o (^p) , (12) 

where that time the scattering amplitude / is not a singular distribution, but a smooth 
function, and the asymptotic behavior can be treated in uniform sense. 

Under our assumptions over the potential the scattered plane waves create for 
E > a complete system of the eigenfunctions of the uniform in multiplicity continuous 
spectrum of the three particle Schrodinger operator, E > 0. 

Our further plan is following: we construct in explicit form a function Y>i(x, q) and 
hope that the difference ip — tpi has the diverging asymptotic behavior 

e »M|q| / 1 \ 
VKx,q) - Vl(x,q) = g(±,q)-^j£ +0 j » ( 13 ) 

where g is a continuous function of the arguments. 

Constructing ipi we use two criteria: (1) The discrepancy 

Qty 1 ](x,q)=-A4 1 + (v(x 1 )+v(x 2 )+v(x!i))ip 1 -Eih, E = |q| 2 . (14) 

sufficiently quickly vanishes at infinity, (2) The asymptotic representation for ip\ — 
e i<x,q> contains asymptotically only the diverging wave. 
Consider the difference 

C = V-^i- (15) 

It satisfies the equation 

H£-E£ = -Q, H = -A + {v(x 1 )+v(x 2 ) + v(x 3 )). (16) 

Since Q is quickly vanishing one can hope that £ asymptotically behaves as the diverg- 
ing wave 

e(x,q)= 5 (x,q)^ + (^J, (17) 
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with a continuous amplitude g. In other words, £ satisfies the classical radiation con- 
ditions at infinity. 

Further, it is naturally to hope that for £ we can construct an integral equation 
with the same properties as the properties of classical Lippmann-Schwinger equation. 
We can do it developing the ideas of work [3]. However, preliminary, we can try to 
use (fTT)|) - (fT7|) for the numerical computation of £ and, consequently, of ip. For the 
numerical computations we can replace (fTTl) by approximate boundary condition 



where R is sufficiently large. 

The following construction of ip\ will consist of two steps. At first, we construct 
for ip\ so called ray approximation tpR. Its discrepancy has some singularities. After a 
natural modification motivated by some classical diffraction problems the discrepancy 
will become a smooth function. 

2.4 Ray approximation 

Consider six vectors <rq. These vectors, more precisely, spanned by them rays, separate 
six sectors that we denote Kp. The indices of the notation coincide with the indices 
of the vector ±1 3 - that belongs to the sector Kj . 

Now we can give explicit expressions for the ray approximation in different sectors 
Kf. 

Sector Kf: ipn = ipf, 

^(x,q) = Xi(x,q)s 2 s 3 . 

We use here the following notations: Sj = s(kj), rj = r(kj). 
Sector K^: ipR = ip^ , 

^(x,q) = X3(x,q)sis 2 . 

Sector : ipn = ip£, 

i>t ( x > °l) = X2(x, q)si + X2(x, T 3 q)s 2 r 3 . 
Sector K^' ipR = ip^i 

?/V(x,q) = X2(x,q)s 3 + X 2 (x,riq)s 2 ri. 
Sector tpR = ip^, 

^r( x >q) = Xi(x,q) +Xi(x,r 2 q)r 2 si + Xi(x, r 3 Tiq)r 2 ri + xi(x, r 3 q)r 3 . 
Sector K^: ipR = rp^, 

^3~(x,q) =X3(x,q) +X3(x,r 2 q)r 2 s 3 + X3(x,Tir 3 q)r 2 r 3 + X3(x,riq)ri. 
The total field ipR is defined by the formula 

TpR = Ot^Pt + W + ^2 + + W + W + «W- 

The notation 6^ is used here for the characteristic function of the corresponding 




for \x\ = R 



(18) 



sector K , 
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In this formula the value of the field i/jr on the boundaries of the sectors is not 
defined. In [27J it was shown that on all boundary rays except two, directed along the 
vectors 

Q23 = t 2 t 3 ci, q 21 = r 2 riq, 

the field is smooth, and its discrepancy everywhere except the two vectors is equal to 
zero. 



2.5 Diffraction corrections 

The diffraction corrections on rays directed along the vectors q 2 3 and q 2 i, can be 
constructed quite easily. Consider the sector A 2 3i containing q 2 3. Introduce the polar 
coordinates (r = |x|,u;). Let us orient the angle from 1^ to 1[\ Let uj 2 % correspond 
to q 2 3. Introduce four angles < u\ < uj 2 < uj 2 % < uj% < CJ4 < tt/3. Consider the 
open covering of the interval (0, vr/3) by the subintervals (0,u; 2 ), (ui,uJi), (u3,7r/3) and 
introduce a subordinated partition of unit: 

1 = Ci + C2 + Cs- (19) 

Further consider the function 

$(a) = ^ I" e it2 dt. (20) 



Notice that 

3>(a) — > 1, as a — > +00, <!>(«) — * 0, as a — ► —00. (21) 

In more detail: 

$(a) = hA$(a), A$(a) = -?-^r ^—zdt = 0(a~ 3 ), when a -> +00. 

Vtt 2ia y7r J a 2ii^ 

Introduce the function 

$S 23) =*(fli<pi(a^3-w)||q23||x|- <q 23 ,x> I 1 / 2 ), (22) 
$( 23) =$( S i 9 n(a;-^ 2 3)||q 2 3||x|- <q 23 ,x> I 1 / 2 ). (23) 



e i < x .q23>$( 23 ) (24) 



It is known that 

</> = e '< x .<i23>^ 

satisfies the Helmholtz equation — A</> — Ecfi = 0. 

Now we can describe the diffraction corrections to the ray approximation on A 2 3i. 
For that the ray field i/jr = 2 ip 2 + ^1^1 ^ n ^ rie sec tor A 2 3i is replaced by 

^g 3) = fa + C2e i< ^ x> [Ri(^ 3) - 6+) + R 2 (^ - ^r)]. (25) 

Ri = ris 2 r 3 , R 2 = r 3 r 2 si + s 3 r 2 r 1 . 

Notice that the field i/jr on the interval (ui,ui4) contains the discontinuous compo- 
nent 

^j = e i< ^ x> [9+Ri + 9^R 2 ], (26) 

so the sense of the modification in nothing else but a simple replacement of this dis- 
continuous on q 2 3 component by a smooth solution of the Helmholtz equation that 
outside of (lj 2 ,Us) gradually transfers to the original discontinuous component up to 
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a diverging circle wave with a smooth amplitude. Outside of the interval (uj\,uj^) the 
function i/ip coincides with the original ray approximation ipR. 

Analogous constructions can be also considered in the sector A312. It is also worth 
to introduce here the polar coordinates, and again to suppose that the angle u varies 
in the same limits with the same orientation, from I3 to — 12. We again can introduce 
the angles u 2 i, cjj, j = 1, 2, 3, 4 and a cutoff function £ 2 . After that the modified field 
on A312 can be described by the formula 

^g 1] =^R + C 2 e i<q21 ' x> [ J R 2 (^S 21) - 0+) + M®? 1 } " 2 )]. (27) 

Here 

$S 21) = <5>{sign{uo 21 - w )||q 2 i||x|- < q 2 i,x> j 1 / 2 ), (28) 

4 2 2l) = ^(sign(u; - u; 21 )||q 21 ||x|- < q 2 i,x > j 1 / 2 ). (29) 

As a result everywhere on T outside of some circle C ri with the center at and the 
radius r\ there appears a smooth approximate wave field ?/>o : 

^0 = ^ R 0l + V>g 3) #231 + ^ #312 ■ (30) 

Here #231 and #312 are the characteristic functions of the corresponding A-sectors, and 
9 1 is the characteristic function of their complement. Again there are no jumps on the 
boundaries of the A-sectors. 

Consider a circle with the center at the origin. The radius n of this circle is defined 
by the condition that outside of the circle on the rays directed along the vectors crq 
the sum of the pair potentials is equal to zero. Under this condition the field ipo can 
be additionally modified with the help of the cutoff function C(l x |) that is equal to 
for |x| < n and to 1 for |x| > r 2 where T\ < r 2 . 

The final expression for the approximate field is now 

V'l = V>oC- (31) 

2.6 Discrepancy 

We remember that there were proposed two criteria that have to be taken into account 
when constructing the function ip\ . It is sufficiently clear that the second one : (2) The 
difference ip\ — e i<x,q> contains asymptotically (in the weak sense) only the diverging 
circle wave, is fulfilled. It is remained to check the first one: (1) The discrepancy 

Qty 1 ](x,q L ) = -Aip 1 + (v{x 1 )+v{x 2 ) + v(x 3 ))Tp 1 -Ei> 1 , E = |q| 2 . (32) 

sufficiently quickly vanishes at infinity. 

From the previous formulas it follows that outside a certain circle of the radius r\ 
the discrepancy is not equal to zero only on some neighborhoods the rays generated by 
the vectors q23 and q2i- On these neighborhoods the discrepancy vanishes as |x|~ 5 / 2 . 
It follows from this that the relative scattering amplitude g(x, q), see (I13p . must be 
continuous. Here we give for the discrepancy a formula that can be used for the 
numerical computations of ip. 

Consider now the field ip\ on the neighborhoods of q23 and q2i. It is not hard to 
see that the discrepancy of this expression is equal to zero on the sectors where there 
is equal to zero the derivative of the function ( 2 . It means that the discrepancy Q[^o] 
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can differ from zero only on the subintervals (o)\,o)2) and (0)3,0:4). That implies that 
the discrepancy on the sector A231 can be naturally represented as the sum: 

= Qf) + Qf) . (33) 

Similarly, on the sector A312 

Q( 21 ) = Qf 1) +Qf ) . (34) 

All four terms here can be easily computed. The answers are completely analogous. 
In particular, 

Qj 23 ) = R x (-A - £)e i<q23 ' x> (^ 23) - l)Ca = (35) 
= R 1 [e i<t * a * > (& { ™ ) - 1)^C2 " < q23,^ > e i<q23 ' x> A$[ 23) ]. (36) 

where w is a unit vector orthogonal to x and oriented along the direction of increasing 

0). 

Finally, 

QbM = Q? 3) + Q? 3} + Q? 1] + Q? 1] ■ (37) 

It is easy to see that all four components of the discrepancy vanish at infinity like 

|x|- 5 / 2 . 

The previous computations of the discrepancy were given for not small |x| where 
the supports of three potentials are separated. Let us modify now the field ip\ by 
introducing in it the factor £ = C(l x l) that is equal to for |x| < r\, and is equal to 
1 for |x| > r2, < ri < r-i- It is supposed that for |x| > r\ three supports do not 
intersect. The definition of ip\ is given by the formula 

V>i = ^oC- (38) 
The final expression for the discrepancy is given by the formula 

Qm = Q[MC - 



d rV'o^q) 



C'-^oA^lxlJV- (39) 
|x| o|x| a|x| 



d|x 

There is no problem in explicit computation of the derivative T^jVo( x ) ^D- 

3 Numerical computations 

The goal of the computations was to show that the suggested plan is realistic and 
can be practically used for the computations of the scattered plane wave and the 
corresponding amplitude of scattering. The pair-particle potential v(x) and the vector 
q are two parameters of the problem. 

As for v (x) we choose the potential function 

M < I (40) 
otherwise. 

Any specific choice is not crucial, we could take arbitrary even potential (even non- 
necessary continuous) with the compact support. With this potential we computed 
the solution x(x, fc) of one-dimensional Schrodinger equation © and found the corre- 
sponding transition s(k) and reflection r(k) coefficients. 
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Then the solutions x{ x i k) were interpolated to the actual computational domain 
to construct the functions Xj- This interpolation was necessary only on the support 
of the potentials. Outside of the supports the analytic expressions of the functions 
x(x, k) were known after the coefficients s(k) and r(k) were found numerically. 

We took E = 4. For the vector q we used two choices: 1) k\ = 1, p\ = \/3, 2) 
k\ = pi = v2- In the first case the field as a function of x is symmetric with respect 
to the straight line generated by q. It was taken for the control. 

The function ipR was computed directly with the knowledge of %ji the Fresnel 
integral was taken from GSL (Gnu scientific library). 

The functions ^o(x,q) and ^/>i(x, q) were computed with the help of the explicit 
formulas for them. The discrepancy Q was also computed with the help of the explicit 
formulas. The radii n < r2 were taken as n =4, r2 = 14.5. We think that this choice 
reasonably corresponds to the selected value of |q|. 

For the diffraction corrections (and near the origin) the chosen partition of unity 
corresponds to function £(z) = z 3 (W — 15z + Gz 2 ), < z < 1, where z is the variable 
relative to angle u. 

Then we finally considered the boundary problem (|16l- 1 18[) . Of course, it was the 
main part of the numerical program of the work. The problem on the disc is not on 
the spectrum. 

For the computations we used mainly FreeFem++, which is a user friendly language 
dedicated for solving partial differential equations with the finite element method. All 
the necessary steps from mesh creation to solving the linear system can be done within 
the same program in a manner that is not of a black box type. Since we used the finite 
element method, we introduced the corresponding weak formulation of the problem: 
find £ G iJ^O) such that 



/ V£ • Vw + (v {x\) + v(x2) + v (xs) — E)^w dx — I i^/E^w dS 
Jn Jon 

= -[Qwdx VwGH 1 ^) (41) 
Jn 



The finite element discretization of (I4ip was then done in a standard fashion using 
quadratic Lagrange elements on a triangular mesh. The computational domain was 
divided into sub-domains to have the finite element mesh fit better with the support of 
the potential V and the constructed function xo and the discrepancy Q. A relatively 
uniform mesh was introduced with lengths of triangle edges between 0.15 and 0.48. 
With a circular domain of radius 190, the total number of degrees of freedom was 3 
million. We used Matlab's solver for large linear systems. 

The results are represented by the Fig. 2-3, and we think that they are reasonable. 



11 



Figure 2 




real(Q) 



Figure 3 
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A certain problem was the choice of the radius R. To get more precise results it 
would be better to take bigger R, but the bigger R means the harder computations. 
The criterium of the compromise was connected with the integral form of the radiation 
conditions. These conditions are: 
l)the integral 

/ ^(x,q)| 2 (42) 
Js R 

over the circle |x| = R must be bounded for large R; 
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2) the integral 

/ ^|(^-z^(x,q)| 2 (43) 
Js R "\x\ 

must decrease as i?~ 2 . 

Notice that Fig. 4 shows that the first integral here is asymptotically approaching 
a constant at sufficiently large |x|, and the second integral is quite small for such |x|, 
but does not decrease for the present computations with R = 190. 




Figure 4 



It, probably, means that such radius is not completely sufficient for the final com- 
putations. However, the bigger radius would mean the harder computations, so we 
decided at the moment to restrict the radius of the circle by 190. 

We considered also the corrected boundary condition where the next term of asymp- 
totic behavior of £ was also taken into account: 

Nevertheless, the correction did not help to stabilize the calculation in smaller domain, 
as it could be expected. The reason is that the term ^ appeared to be very small 
comparatively with other terms. 

To clarify further the situation with the stabilization of L2 - norm on the boundary 
of the disk we should come back to behavior of £ as a function of the angle at fixed 
radius r = R. 
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Figure 5 



One can easily see from the Figure that the main contribution to the integral 
comes from two special directions on the boundaries of shadow and light in sectors 
A231 and A312, see Fig. 5. Therefore for the stabilization of the whole integrals first of 
all the contributions to them of two indicated sectors must be stabilized. Namely their 
stabilization was not completely reached for considered size of configuration domain 
and requires bigger scale of radii. There are also other computational indications that 
the behavior of the field in these two sectors is responsible for the (non)stabilization 
of the L2 norms. 

On the other hand, we found the right tendency of the behavior of the solution for 
large r what was the aim of the present calculations. 
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